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METHOD AND APPARATUS FOR 
CONTROLLING ABLATION IN REFRACTIVE SURGERY 

T ' ' ■ — - . . . — 

Background of the Invention 

Laser refractive surgery is often used to correct 
5 refractive errors, such as myopia, hyperopia, and 

astigmatism, in a patient's eyes. Refractive errors 
are also called lower-order aberrations of the eye. 
The eye also may have higher-order aberrations caused 
by irregularities in the cornea or crystalline lens. 

10 More recently, laser correction is also applied to 

higher-order aberrations such as coma, spherical 
aberrations, and other aberrations. Ideally, the laser 
refractive surgery should accurately remove the 
refractive error and. not induce additional aberrations. 

15 The refractive errors and aberrations are removed 

through corneal ablation. The laser device follows a 
predetermined ablation pattern designed to correct the 
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refractive errors and aberrations. These patterns are 
referred to as ablation patterns. 

If the laser ablation patterns used in 
photorefractive keratectomy (PRK) and laser in-situ 
5 keratomileusis (LAS IK) procedures to correct the 

refractive error are not properly selected, 
post-procedure aberrations can result that reduce 
quality of vision in scotopic conditions. These 
resulting post-procedure aberrations . are not inherent 

10 in the predetermined ablation patterns used to correct 

the refractive error, but are believed to result from 
the corneal healing response after surgery. Such 
post-procedure or post-operative aberrations are 
referred to herein as "induced aberrations." 

15 Known laser ablation algorithms for vision 

correction have assumed an ideal target spherical 
surface. This is known in the art to be the Munnerlyn 
approach. The pre-operative corneal shape is 
approximated by the closest spherocylindrical shape. 

20 The ablation pattern is the difference between the 

assumed ideal target shape and the pre-operative 
spherocylindrical shape. The Munnerlyn approach 
specifies that the ablation pattern is limited to an 
area referred to as the "optical zone/' or "OZ". where 
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full correction is intended. Surrounding the optical 
zone is an area referred to as the "transition zone" 01 
"TZ " The transition zone is added to avoid an abrupt 
transition in the shape of the cornea at the edge of 
5 the optical zone. Prior known transition zone designs 

have provided continuity of ablation depth but do not 
guarantee continuity in slope and curvature. 

Known formulations of ablation patterns are also 
algebraically complex in that they result from 

10 differences between spherical surfaces. However, 

spherical surface formulations contain spherical 
aberrations. In addition, spherical surfaces are not 
ideal refractive surfaces. 

Ablation patterns have been proposed that 

15 compensate for induced spherical aberrations using an 

empirical approach. However, these proposed ablation 
patterns only correct for spherical myopia and are not 
useful to correct myopic astigmatism, hyperopia, 
hyperopic astigmatism, or mixed astigmatism. 

20 Elliptical transition zones for astigmatic 

ablations have been proposed which increase the 
transitional zone for ablations with higher slope 
transition. This approach limits the slope of the 
transition zone. 
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Wavef ront-guided laser vision correction 
procedures have been also used. This procedure 
measures the aberration of the eye and attempts to 
correct them. Theoretically, this should decrease the 
5 aberration of the eye after surgery. However, it is 

believed that this procedure results in secondary 
aberrations greater than the aberration corrected. 

Summary of the Invention 

The present invention relates to a method and 
10 apparatus for controlling ablation in refractive 

surgery using laser ablation patterns to correct 
refractive errors of the eye such as nearsightedness 
(myopia) , farsightedness (hyperopia) , astigmatism, and 
higher order aberrations of the eye. In accordance 
15 with the present invention, post-procedure (also 

referred to as "post-ablation" ) surface smoothing is 
predicted and the ablation pattern modified to based on 
the prediction to correct for induced aberrations. 

A method is provided for designing laser ablation 
20 patterns that corrects refractive errors and 

compensates for aberrations from post-ablation corneal 
surface smoothing function. A mathematical model is 
used to characterize this smoothing function and 
compensate for post-operative changes to prevent 
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induceci optical abe rrations . The method also controls 
positioning of laser pulses to optimize the achievement 
of the intended ablation pattern. 

In accordance with one exemplary embodiment of the 
5 present invention, an apparatus is provided for 

correcting optical aberrations, the apparatus 
comprising a laser for performing ablation to modify 
the optical aberrations of the patient's eye, and 
control means for controlling ablation by the laser, in 

10 accordance with a control algorithm that corrects for 

the optical aberrations and further corrects for 
predicted induced aberrations. 

In accordance with a preferred embodiment, the 
control means includes means for storing a mathematical 

15 model that predicts the induced aberrations resulting 

from surface smoothing, the control means controlling 
the laser ablation to correct existing refractive error 
and aberration plus induced aberration predicted by the 
mathematical model. The adjustment for induced 

20 aberration is derived from the mathematical model 

preferably using a convolution algorithm to predict 
corneal smoothing. The control algorithm preferably 
provides a transition zone around an optical zone in 
accordance with a cubic spline function. The control 
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algorithm preferably corrects myopia and hyperopia by 
targeting a parabolic-shaped corneal height change. 
The control means further includes means for storing a 
rain-drop algorithm to control pulse placement of the 
5 laser to control ablation. The control means 

preferably includes a sorting algorithm to order a 
laser pulse sequence to minimize the spatial overlap 
between consecutive laser pulses. 

In accordance with another exemplary embodiment of 

10 the present invention, a method is provided for 

modifying refraction of a patient's eye comprising the 
steps of providing a vision correction laser to modify 
refraction of the patient's eye, and controlling 
ablation by said laser in accordance with a control 

15 algorithm that corrects for existing measured 

aberrations and anticipates and corrects for induced 
aberrations . 

In accordance with a preferred embodiment of the 
method, the step of controlling includes determining a 
20 correction map in response to the preoperative 

measurements, storing a mathematical model that 
predicts corneal surface smoothing responsive to laser 
ablation, and determining an ablation map by adjusting 
the correction map in response to the mathematical ' 
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model. Preferably, the step of determining the 
ablation map further includes generating ablation 
patterns using an iterative deconvolution algorithm 
that compensates for a predicted surface smoothing 
5 response to said ablation. Preferably, the step of 

controlling further includes the step using a rain-drop 
algorithm to control pulse placement of said laser to 
control ablation. 

In accordance with yet another exemplary 

10 embodiment of the present invention, a computer program 

product is provided operative in a laser surgical 
device for correcting optical aberrations, the computer 
program product comprising a correction mapping stage 
that determines a correction map based on initial 

15 optical measurements, an ablation mapping stage that 

uses a constrained iterative deconvolution algorithm to 
compute an ablation map from the correction map, and a 
pulse control stage using a second iterative algorithm 
to optimize placement of pulses to avoid spatial 

20 overlap. 

In accordance with yet another exemplary 
embodiment of the present invention, a method for 
correcting optical aberrations comprising the steps of 
providing a computer controlled ablation device for 
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ablating a patient's eye, performing pre-operative 
optical measurements to establish a corrective 
prescription, establishing a mathematical model to 
predict surface smoothing, establishing an ablation map 
5 based on the corrective prescription and the 

mathematical model, and ablating the patient's eye in 
accordance with the ablation map. 

In accordance with the present invention, ablation 
patterns or designs compensate for the predicted 

10 effects of surface smoothing and thereby improve the 

accuracy of refractive correction and reduce 
undesirable secondary or induced aberrations due to 
corneal healing. The ablation patterns in accordance 
with the present invention can be used to correct the 

15 full range of refractive errors such as myopia 

(nearsightedness) , hyperopia (farsightedness) , 
astigmatism, and higher order aberrations. The method 
and apparatus of the present invention provides: 

1. a mathematical model to predict corneal 
surface smoothing response to laser ablation; 

2. a deconvolution algorithm to generate 
ablation patterns that pre-compensate for the 
surface smoothing response; 

3. a healing-adjusted ablation pattern including 
a transition zone around the optical zone for 
providing continuity in ablation depth and 
slope; 
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ablation patterns to correct hyperopia, 
myopia, and astigmatism based on parabolic 
surfaces in accordance with Zernike 
polynomials rather than spherocylindrical 
surfaces (the ideal refractive surface is 
parabolic) ; 

ablation patterns to correct higher-order 
aberration along with the lower-order 
refractive errors ; 

an iterative algorithm to generate a sequence 
of laser pulse placement that produces a 
result very closely matching the desired 
ablation map (The sequence has a fractal 
property such that if it is interrupted at 
any point, the achieved ablation will still 
approximate the shape of the complete 
pattern) ; and 

a sorting algorithm to order the pulse 
sequence to minimize spatial overlap between 
consecutive or temporally nearby laser 
pulses . 



In accordance with the present invention, a 
control algorithm produces a more accurate outcome and 
reduces induced aberrations for the full range of 
5 correction of myopia, hyperopia, astigmatism, and 

higher-order aberration. 

Brief Description of the Drawings 

Fig. 1 depicts an optical zone ("OZ") surrounded 
by a transition zone ("TZ"); 
10 Fig. 2A depicts correction and ablation profiles 

for myopia; 
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Fig. 2B depicts correction and ablation profiles 
for hyperopia; 

Fig. 3A depicts correction and ablation profiles 
for myopic astigmatism with flat meridian (-1.50D 
5 sphere, -1.00D cylinder x 180 ° ) ; 

Fig. 3B depicts correction and ablation profiles 
for myopic astigmatism with steep meridian (-1.50D 
sphere, -1.00D cylinder x 180°); 

Fig. 4A depicts ablation maps for myopia (-1.00D, 
10 OZ=6.0mm, TZ=0 . 6mm) ; 

Fig. 4B depicts an ablation map for hyperopia 
(+1.00D, OZ=5.5mm, TZ=1.8mm); 

Fig. 4C depicts an ablation map for myopic 
astigmatism (-1.50D sphere, -1.00D cylinder x 180°, 
15 OZ=6.0mm); 

Fig. 4D depicts an ablation map for hyperopic 
astigmatism (+O.O0D sphere, +1.00D cylinder x 180°, 
OZ=5 . 5mm, TZ=1 . 8mm) ; 

Fig. 5A depicts an ablation pulse map for 1 D of 
2 0 myopia; 

Fig. 5B depicts an ablation raw remainder map (the 
difference between achieved and target ablation maps) ; 
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Fig. 5C depicts an ablation residual map after 
surface smoothing (difference between achieved and 
target correction maps) ; 

Fig. 6 is a functional block diagram of a laser 
5 refractive surgical device in accordance with the 

present invention; 

Fig. 7 is a schematic view of an eye; and 

Figs. 8-11 are flow diagrams of the ablation 
process in accordance with the present invention; 
10 Fig. 12 depicts an ablation profile for 1 D of 

hyperopic correction with 5.5 mm diameter optical zone 
and a 1.8 mm wide transition zone (the achieved surface 
height change as predicted by the surface smoothing 
model is shown in dotted line and the difference 
15 between the dotted and solid lines is attributed to 

changes in the epithelial thickness) ; 

Fig. 13 depicts an ablation profile for 1 D of 
myopic correction with 6 . 0 mm diameter optical zone; 

Fig. 14A depicts a simulation of the ablation 
20 pattern for minus-cylinder correction on a EC-5000 

with 5 . 5 mm diameter optical zone and 7.0 mm diameter 
transition zone with flat meridian; 

Fig. 14B depicts a simulation of the ablation 
pattern for minus-cylinder correction on a EC-5000 
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with 5.5 mm diameter optical zone and 7.0 mm diameter 
transition zone with steep meridian; 

Fig. 15A depicts clinical results of spherical 
hyperopia correction with slope = -0.708 ± 0.045 
5 (mean ± standard error) with intercept set to zero on a 

LADARVi s i on System (SIRC SE = surgically-induce 
refractive change in spherical equivalent and Laser 
SE = laser ablation for spherical equivalent) ; 

Fig. 15B depicts clinical results of spherical 
10 myopia correction with slope = -0.968 ± 0.014. (mean ± 

standard error) with intercept set to zero on a 
LADARVision System (SIRC SE = surgically-induce 
refractive change in spherical equivalent and Laser SE 
= laser ablation for spherical equivalent); 
15 Fig. 16A depicts correction/ablation ratios as 

functions of the smoothing constant s for simulated 1 D 
ablations, of hyperopia; 

Fig. 16B depicts correction/ablation ratios as 
functions of the smoothing constant s for simulated 1 D 
20 ablations of myopia; 

Fig. 17A depicts an induced spherical aberration 
(coefficient for Zernike series term Z 4 °) as functions 
of the smoothing constant s for simulated 1 D ablations 
of hyperopia; 
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Fig. 17B depicts an induced spherical aberration 
(coefficient for 2 4 °) as functions of the smoothing 
constant s for simulated 1 D ablations of myopia; 

Fig. 18A depicts model simulation of a 1 D 
5 against-the-rule minus cylinder ablation on the EC-5000 

(the flat meridian is vertical and the steep meridian 
is horizontal) showing the ablation map; 

Fig. 18B depicts model simulation of a 1 D 
against-the-rule minus cylinder ablation on the EC-5000 
10 (the flat meridian is vertical and the steep meridian 



is horizontal) showing the astigmatism 



correction/ablation ratio as a function of the 



smoothing constant s; and 



Fig. 18C depicts model simulation of a 1 D 



15 



against-the-rule minus cylinder ablation on the EC-5000 



(the flat meridian is vertical and the steep meridian 



is horizontal) showing surgically-induced epithelial 



thickness change with s = 0.5 mm. 



Detailed Description of Invention 



20 



Ablation designs, in accordance with the present 



invention, start with a target corneal surface height 



change needed to correct preoperative refractive error 



or aberration. 



This is referred to as the correction 



map Ah. The correction map is specified within a 
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central optical zone ( VN OZ") . The OZ preferably centers 
on the line of sight and matches the maximum size of 
the pupil. In accordance with the present invention, 
the correction map for the correction of myopia, 
5 hyperopia, and astigmatism are parabolic. This reduces 

the induction of aberrations compared with spherical 
and cylindrical corrections. 

Referring to Fig. 1, the OZ is surrounded by a 
transition zone ( XN TZ") to produce a smooth, continuous 

10 blend in ablation depth, slope-, and curvature with the 

surrounding cornea. The OZ and TZ together constitute 
the entire ablation zone ("AZ") . In accordance with 
the present invention, a constrained iterative 
deconvolution algorithm is used to compute the ablation 

15 map for the AZ, which, after the expected postoperative 

surface smoothing effects occur, produces the desired 
correction in the OZ . 

The collection of laser pulses needed to produce 
the ablation- map is then computed in a second iterative 

20 algorithm that optimizes the placement of laser pulses. 

The pulse sequence is then sorted to avoid spatial 
overlap between consecutive or temporally nearby 
pulses . 
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In general, the refractive error correction for an 
eye must correct for defocus (hyperopia or myopia) , 
astigmatism, and higher order aberrations. The 
correction maps, in accordance with the present 
5 invention, for all three are described below. Further, 

in accordance with the present invention, these 
correction maps can be combined to customize the 
ablation for an individual eye. The ablation map and 
pulse sequence are then generated using the iterative 
10 algorithms. 

Correction Map Generation 

Correction Map for Defocus 

Defocus of the eye is caused by the focusing power 
of the eye being too strong (myopia or nearsightedness) 

15 or too weak (hyperopia or farsightedness) . To correct 

defocus without adding aberration, the ideal corrective 
lens should be parabolic rather than spherical. 
Spherical surfaces are associated with spherical 
aberration. Spherical lens have been historically used 

20 in spectacle correction, not for optical superiority, 

but because traditional lens grinding methods produce 
spherical surfaces much more readily than aspheric 
surfaces . 
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With laser refractive surgery, aspheric ablation 
can be performed just as easily as spherical patterns. 
The. ablation designs in accordance with the present 
invention are based on a more optically correct 
5 parabolic shape. 

Optical surfaces with parabolic shapes are based 
on the second order of the Zernike circle polynomial 
series. The Zernike series is a standard for analyzing 
optical aberration and optical surfaces. In Zernike 

10 series, the surface is defined within a circle of unit 

radius. Each Zernike series term Z n m has an order n and 
an angular frequency m. Defocus (myopia and hyperopia) 
is described by Z 2 ° and astigmatism is described by Z 2 ±2 . 
The preferred correction maps of the present invention 

15 for defocus (expressed by Equation 1 below) are defined 

by the defocus term Z 2 °(p,6) (expressed by. Equation 2 
below) . 

Defocus error of the eye is conventionally 
measured by "manifest refraction." Manifest refraction 
20 is performed by trial correction with spherical and 

cylindrical lenses. The combination of lenses that 
receive the best subjective rating is the manifest 
refraction, which can be directly translated to a 
prescription for spectacle lenses. Spectacle 
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refraction is generally written in notation such as 
"sphere cylinder * axis." "Sphere" is the power of the 
spherical lens, "cylinder" is the power of the • 
cylindrical lens, and "axis" is the orientation of the 
5 cylindrical lens. Defocus is measured by the spherical 

equivalent ("SE") , which is equal to sphere + 0.5 * 
cylinder. A formula (expressed by Equation 3 below) is 
provided for converting the SE to the Zernike 
coefficient used in Equation 1 below. The corneal 
10 index is substituted for the keratometric index. 



Defocus can be measured objectively by a wavefront 
sensor or an autoref ractor . With a wavefront deviation 
measurement, Zernike decomposition is usually 
determined by analysis software. The defocus 



15 



coefficient w 2 ,o from the wavefront decomposition 



specifies the required defocus correction ( expressed 



by Equation 4 below) . 



^M)=AA 2 ,oZ 2 0 (A^) 



Equation 1 



20 




Equation 2 



A* 2i0 =-/)„*.7[W3(n-l)] 



Equation 3 
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A/Z 2,0 = "W(*-l) 

Equation 4 

where 

Ah De focus (r, 0) = the correction map (target corneal 
5 surface height change) in pm, 

r = radius from the center of the optical zone in 
mm, perpendicular to the line of sight, 

9 = meridian angle which is zero along +x (to the 
right when facing the front of the eye) and increases 
10 counterclockwise, this follows the right-handed 

coordinate convention with the height dimension z 
positive outward along the line of sight, 
p = r/Roz is the normalized radius, 
Roz = radius of optical zone in mm, 
15 Ah 2 ,o = coefficient for the Z 2 ° term of the 

correction, map, 

Z2 0 (p, 9) = Zernike polynomial for defocus, 
Dse = spherical equivalent correction in diopter 
(myopia negative, hyperopia positive) specified at the 
20 corneal plane, 

n = corneal refractive index, a value of 1.377 is 
used, and 

w 2 , 0 = coefficient for the Z 2 ° term of wavefront 
deviation. 
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For myopia correction, the diameter of the OZ is 
preferably 6 . 0 mm but adjustable between 5.5 mm 
and 6.5 mm. Hyperopia correction requires a larger TZ 
and therefore the OZ has to be relatively smaller to 
5 fit inside a reasonable AZ . For hyperopia correction, 

the diameter of the OZ is preferably 5.5 mm but 
adjustable between 5 . 0 mm and 6.0 mm. A smaller OZ may 
be needed to fit under a smaller LAS IK (laser in-situ 
keratomileusis) flap or to reduce ablation depth in 

10 cases of large corrections and thin corneas. A larger 

OZ may be needed in eyes with larger pupils. 

The correction profiles for 1 diopter ("D") of 
myopia and hyperopia are shown on Fig. 2A and 2B, 
respectively. The corresponding ablations maps are 

15 shown on Fig. 4A and 4B. 

Astigmatism Correction Map 

The astigmatism correction (expressed by 
Equation 5 below) is preferably defined with the 
Zernike terms Z 2 2 (p,.A) and Z 2 ~ 2 (p,A), Equation 6 below, 
20 which represent cardinal and oblique astigmatism, 

respectively. The diameter of the OZ is preferably 
matched to that used for the correction of defocus. 

When astigmatism is measured by manifest 
refraction, Equation 7 below is used to obtain the 
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Zernike coefficients for the correction map. 
Astigmatism can also be measured by a wavefront sensor 
which can be converted to the correction coefficients 
with Equation 8 below. Finally, if astigmatism is 
5 measured by corneal topography, it can be directly 

translated into the correction coefficients (Equation 
below) . 

Ah Ast ig(r,e) = Ah 2,2 Z 2 2 (p,9) + Ah 2 ,-2 Z 2 ~ 2 (p,9) 

Equation 5 

10 Z 2 \p y 9) = 46p 2 cos2e 

Equation 6A 

Z 2 -\p,0) = y/6p 2 sin 20 

Equation 6B 
= D asus*oz 2 cos(2 x axis)l[A46{n - 1)] 
15 Equation 7A 

A V 2 = D asu Z Kz sin(2 x axis)/[4yf6(n - 1)] 

Equation 7B 

Ah 2,2 = w 2 , 2 /(n -1) 

Equation 8A 

20 Ah 2 ,-2 = w 2 ,_ 2 /(n -1) 

Equation 8B 

Ah 2,2 = ~h 2 ,2 

Equation 9A 
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Ah 2,-2 = ~h 2/ -2 

Equation 9B 

where 

Ah Astig (r f 0) = the correction map (target corneal 
5 surface height change) in ym, 

Z 2 2 (p,0) = Zernike polynomial for cardinal 
astigmatism, 

Z.2~ 2 (p,6) = Zernike polynomial for oblique 
astigmatism, 

10 Ah 2 , ±2 = coefficients for the Z 2 ±2 terms of 

correction map, 

^Astig — astigmatism magnitude in diopters from 
refraction at the corneal plane, 
axis = astigmatism axis, 
15 w 2/±2 = coefficients for the Z 2 ±2 terms of wavefront 

deviation, and 

h 2 ,±2 = coefficients for the Z 2 ±2 terms of corneal 
topographic height. 

2 —2 

Z 2 and Z 2 contain no defocus power and correspond 
20 to Jackson cross cylinders (except that the Zernike 

terms have parabolic profile and the cylinder has 
circular profile, which is slightly different) . 



WO 03/075778 PCT/USO J/06343 

-22- 

Higher Order Aberration Correction Map 

Higher order aberrations are those that are of 
more complex shape than defocus and astigmatism. 
Higher order aberration is described by Zernike terms 
5 of order greater than 2. For example, coma is 

primarily described by Z 2 ±l and spherical aberration is 
primarily described by Z 4 °. Higher order aberration can 
be measured by either corneal topography or. a wavefront 
sensor. They are separated from lower order 
10 aberrations using Zernike series decomposition 

(expressed by Equation 10 below) . 
h H o(p,e) = h(p,6) - h 2 , 2 Z 2 2 (p,6)' - h 2 ,- 2 Z 2 - 2 (p, 9) - 

h 2 ,oZ2°(P,e)- hi^z^tp.e) - h^zr 1 (p,e) - h 0 ,oZo°(p,e> 

Equation 10A 

15 

w H0 (P,e) = w(p,6) - w 2f2 Z 2 2 (p f 0) - w 2 ,- 2 Z 2 ~ 2 (p,9) " 
w 2 , 0 Z 2 °(p,e)- w^xZ^Cp^e) - w^^zrMprS) - w 0 , 0 Z 0 °(p,e) 

Equation 10B 

20 where 

h(p,0) is the corneal topographic height, 
subscript H o denotes sum of higher-order terms, 
p = r/R oz , the normalized radius, 
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^n f m are the coefficients from Zernike 
decomposition of h(p,9), 

w(p,Q) iis the wavefront height, and 
w n , m are the coefficients from Zernike 
5 decomposition of w(p,9). 

The correction map &h H0 (r,e) for higher order 
aberration is described by Equation 11 (topography- 
guided) or Equation 12 below ( wavef ront-guided) . 

Ah HO (r, 6) = -h HO (r, 9) 
10 Equation 11 

Ah HO (r,e) - w HO (r,9) / (n-1) 

Equation 12 

Combining Correction Maps 

The correction maps for defocus and astigmatism 
15 are combined to produce the target refractive 

correction. If higher order aberration has been 
measured by a wavefront sensor or corneal topography 
system, it can also be added to the target correction 
map at this step. All of the component maps must have 
20 the same OZ diameter. The combination is described in 

Equation 13A below. The maximum point of the 
correction map is adjusted to zero in Equation 13B 
below. 

Ah temp (r, 9) = Ahoefocus (r/9) + Ah Astig (r,e) + Ah HO (r,9) 
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for r < R oz 

Equation 13A 
Ahtarget (r,0) = Ah temp (r / G) - max{ Ah temp (r,e) } 
for r < Roz 

5 Equation 13B 

where 

Ah t emp(r,8) is the combined target correction map 
before ad j ustment , 

Roz is the optical zone radius in global 
10 coordinates, 

Ahtarget (r, 9) is the final combined target 
correction map, and 

max{} returns the maximum height of the enclosed 
surface . 

15 Ablation Map Generation with Iterative 

Deconvolution Epithelial Smoothing Model 

The corneal surface change specified by * 

Ah target (x, y) will be different from the ablation map 

a(x,y) because of a corneal surface smoothing response. 

20 One mechanism of smoothing is the migration of surface 

epithelial cells away from more convex areas (islands) 
and into less convex areas (divots) . In accordance 
with the present invention, this smoothing response can 
be represented by a mathematical model that predicts 

25 this smoothing response. Details regarding 



10 
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construction of the mathematical model to predict the 
smoothing response in accordance with- the present 
invention are described below in detail under the 
section entitled "Construction of Mathematical Model of 
5 Corneal Surface Smoothing after Laser Refractive 

Surgery . " 

According to the mathematical model of the present 
invention, Ah'(x,y) relates to a'(x,y) by a convolution 
operation (expressed by Equation 14A below) with 
smoothing function f'{x' ,y') (expressed by Equation 14C 
below) . 

A*Vj>aVj , )®/'(x , ,y) 

Equation 14A 
AH'(co x ', CD/) =A'(g> x ', co/) F'(cd/, CD/) 

Equation 14B 

F'(CD X ', CD/) = 1/[1 + S 2 (CD X ' 2 + CD/ 2 )] 

Equation 14C 

where 

Ah' (x',y') is the corneal surface height change in 
20 local coordinates, 

a'(x',y') is the ablation map in local 
coordinates, 

f'(x',y') is the impulse response function of the 
corneal surface (convolution with f (x',y') describes 



15 
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corneal surface smoothing, f'(x',y') is the inverse 
Fourier transform of F' (co x ' , co y ' ) ) , 

AH' (co x ' , cd y ' ) is the Fourier transform of 
Ah' (x' ,y' ) , 

5 u x ' and co/ are the respective spatial frequencies 

. for x' and y' in radians/length, 

A' (co x ' , co y ' ) is the Fourier transform of a' (x' , y' ) , 
F' (co x ' , co y ' ) is frequency response of the corneal 
surface (it's the Fourier Transform of f (x',y')), and 
0 s is the smoothing constant that characterizes the 

epithelial smoothing model (1/s is the cutoff frequency 
(radian/length) of the low-pass filter (Equation 14C) ) . 

Equations 14A and 14B are equivalent expressions 
of the corneal surface smoothing model in the spatial 
5 and frequency domains, respectively. The spatial and 

frequency domains are related by the Fourier transform. 
The frequency response of the corneal surface 
(Equation 14C) has the form of two first-order 
Butterworth low-pass filters in series. The filter is 
D characterized by a smoothing constant s. Clinical data 

from laser in-situ keratomileusis (LAS IK) has been used 
to estimate the value of s. The best-fit values for 
corrections of myopia, myopic astigmatism, and 
hyperopia have been found to range between 0.32 
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and 0.63 mm. . The average was 0.5 mm. According to 
these results, .s = 0.5 mm is preferably used to 
generate the ablation map in the deconvolution 
. algorithm described below. However, the control 
5 ..algorithm of the present invention is not restricted to 
a specific value of smoothing constant. 

Conversion Between Local and Global Coordinate Systems 

Equations 14A-C are expressed in local coordinate 
system x' , y' which are tangential to the corneal 

10 surface. In the local coordinate system, z' is defined 

as being perpendicular to the local corneal surface. 
The local coordinate system has the perspective of a 
small epithelial cell on the corneal surface and is 
important to the development of the epithelial 

15 smoothing theory. The corneal height and ablation 

depth are preferably specified in fixed global 
coordinate system x, y, z. In the global coordinate 
system, z is defined as the line of sight. The global 
coordinate system has the perspective of the laser 

20 scanning system looking down at the cornea. The 

coordinate systems are identical at the corneal apex 
but the deviation increases as the radial distance 
approaches the corneal radius of curvature. 
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The following approximate coordinate 
conversions are used to apply the smoothing theory 
(Equations 14A-C) to ablation design. The equations 
treat the corneal surface as a perturbation from a 
5 sphere. 

a'(r',0) = a(r,9) cos [arcsin (r/R c ) ] 

Equation 15A 
a (r, 9) = a' (r' , 9) /cos (r # /R c ) 

Equation 15B 

10 r' = R c arcsin(r/R c ) 

Equation 15C 

r = R c sin(r' /R c ) 

Equation 15D 
Ah'(r',9) = Ah(r,9) cos [arcsin (r/ R c ) ] 
15 Equation 15E 

Ah(r,9) = Ah'(r\0)/cos(r'/R c ) 

Equation 15F 

where 

r' is the local radial coordinate equivalent to 
20 the arc length from the corneal apex (the apex is 

defined as the intersection between the line of sight 
and the anterior corneal surface) ; 

r is the global radial coordinate, the radial 
distance measured perpendicularly from the line of 
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sight (the line of sight is defined as the line passing 
through the center of the eye's entrance and exit 
pupils connecting the object of regard to the foveola) ; 

R c is the radius of curvature of the corneal 
5 best-fit sphere (a value of 7 . 6 mm is used in 

simulations ) ; 

global height guantities h, Ah, and a are measured 
along the z axis, which is positive outward along the 
line of sight; and 
10 6 is the meridian angle or the azimuth, which is 

zero on the +x axis and increases counterclockwise. 

Conversion Between Polar and 
Cartesian Coordinate Systems 

The coordinate transforms in Equations 15A-F. are 

15 performed in polar coordinates. To convert to and from 

the Cartesian coordinate system, the following 

equations are used. 

x = r cos 0 



20 y = r sin 8 



r = y/{* 2 + y 2 ) 



6 = arctan (y/x) 



Equation 16A 



Equation 16B 



Equation 16C 
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10 



20 



Equation 16D 



and 



x' = r' cos 9 



y' = r' sin 9 



Equation 17A 



Equation 17B 



Equation 17C 



G = Cretan (y' /x' ) 

Equation. 17D 
where - 

x is perpendicular to the line of sight and 
positive to the right of a viewer facing the eye; and 
y is perpendicular to the line of sight and 
15 positive upward (with face upright) . 

Deconvolution To Compute Expected Smoothing 

The ablation pattern that compensates for the 
expected smoothing change can be computed by 
deconvolution. Deconvolution is the reverse of the 
convolution operation shown in Equation 14A. It can be 
written in the frequency domain as: 

A'(g> x ', co/) - AH' (q/, Q y ' ) / F'(co x ', co y ' ) 

Equation 18 
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Equation 18 is generally not solved directly 
because F' approaches zero at high spatial frequencies. 
Division by zero is not defined. Division by a very 
small number is computationally susceptible to noise or 
5 rounding errors. Stable solution is usually obtained 

by constrained iterative deconvolution where the 
smoothing change described by Equation 14A. is 
subtracted from the ablation pattern until the ideal 
ablation map converges. 

10 Transition Zone Width Determination 

The surface smoothing model, in accordance with 
the present invention, indicates that any sudden change 
in depth or slope on the corneal surface would be 
smoothed over. To obtain the desired surface change in 
15 the optical zone, it must be. surrounded with a gradual 

transition. The transition zone occupies the area 
between the borders of the ablation zone and the 
optical zone. 

a'(r',6) = spline' (r' ,9) for R AZ ' (9) > r' > R 02 
20 Equation 19A 

W TZ ' (9) = R A2 ' (9) - R 0 z' 

Equation 19B 

W TZ (9) = RazO) - R 0 z 

Equation 19C 
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where 

spline' (r',9) is the spline function for the TZ, 
Roz' is the optical zone radius in local polar 
coordinates , 

5 R A z' (9) is the ablation zone radius in local polar 

coordinates, 

W TZ ' (9) is the transition zone width in local pola 
coordinates, 

Roz is the optical zone radius in global polar 
10 coordinates, and 

RazO) is the ablation zone radius in global polar 
coordinates, and 

W TZ (9) is the transition zone width in global pola 
coordinates . 

15 The width W TZ ' of the transition zone (TZ) is 

preferably set so that the transition spline profile 
has component frequencies mostly below that of the 
cutoff radial frequency 1/s radian/mm. Hyperopic 
correction profiles need a wider transition zone 

20 because the transition contains more phases (Fig. 2B) 

compared to a myopic correction (Fig. 2A) . The 
transition zone for the myopia correction profile 
(Fig. 2A) goes through roughly 1/4 cycle of a sinusoid 
(n/2 radian) . The transition zone for the hyperopia 
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correction profile (Fig. 2B) goes through roughly 3/4 
cycle of a sinusoid (3n/2 radian) . Therefore, W T2 ' 
should be approximately 1/2 n s for myopia and 1.5 n s 
for hyperopia correction. For s = 0.5 mm, 
5 W TZ ' = 1/2 n s translates to W TZ ' = 0 . 8 mm for myopia 

correction. Converted into global coordinates, W TZ = 
0.6 mm. Proportionally, the preferred W TZ = 1.8 mm for 
hyperopia corrections . 

The TZ for myopic astigmatism is more complicated. 

10 On the steep meridian, the TZ is identical to that for 

myopia and W TZ = 0.6 mm is used. On the flat meridian, 
the TZ profile gains more phase with increasing ratio 
of astigmatism to myopia. Therefore, the preferred 
shape of the AZ is elliptical and the width of the TZ 

15 depends on the meridian. The widths on the flat and 

steep meridians determine the diameters of the 
elliptical AZ on these cardinal meridians. The 
diameters for the AZ and the TZ widths can then be 
calculated on all meridians using the equations for an 

20 ellipse. 

For abs(D SE ) > 0.5 D Ast ig and D SE < 0 (myopic 
astigmatism) 

W T2 _ steep = 0.6 mm, 

Equation 20A1 
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W T z_fiat = W TZ _ steep + min(l. 8, R 0 z {sqrt[(-D SE +0.5 D Astig )/- 

D SE - 0.5 D Astig ) ] - 1} ) 

Equation 20A2 

where 

5 min() returns the minimum member of the input set, 

abs() is the absolute value function, 
sqrt() is square root function, 

D S e = spherical equivalent correction in diopter 
(myopia negative, hyperopia positive) specified at the 
10 corneal plane, and 

D A stig = astigmatism magnitude in diopters 
(positive) from refraction at the corneal plane. 

For mixed astigmatism, the flat meridian has a 
profile similar to that for hyperopia and a 1.8 mm TZ 
15 width is preferred. The transition profile on the 

steep meridian depends on the SE. Therefore, the 
preferred AZ shape is again elliptical. In a preferred 
system, the steep meridian's TZ width transitions 
from 0.6 mm to 1.8 mm, depending on the ratio of SE to 
20 astigmatism, according to Equations 20B1-2. 

For abs(D SE ) < 0.5 D As tig (mixed astigmatism) 

W T z_steep = 1.2 + (1.2 D SE / D Astig ) 

Equation 20B1 

^Tz_f iat — 1-8 mm 
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Equation 20B2 
For mixed astigmatism, all meridians have profiles 
similar to that for hyperopia and a 1 . 8 mm TZ width is 
preferred. 

5 For abs(D SE ) > 0.5 D Astig and D SE > 0 (hyperopic 

astigmatism) 

W TZ _ steep = 1.8 mm 

Equation 20C1 

WTz__fiat — 1-8 mm 

10 Equation 20C2 

In most cases, the TZ width does not need to be 
adjusted for higher order aberration because it is 
generally much smaller in magnitude relative to the 
coexisting defocus and astigmatism. However, if the 

15 higher order aberration is larger than the defocus and 

astigmatism, then the larger 1 . 8 mm TZ width is 
preferred all around. 

Cubic Spline In The Ti:ansit:ion Zone 

A cubic spline curve in the TZ that provides a 
20 continuous transition in depth and radial slope is 

preferred. Discontinuity in depth and slope provokes 
large healing responses (epithelial hyperoplasia, 
subepithelial haze) and cause tear film instability. 
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The cubic spline function is a 3 th order polynomial 
in r' (Equation 21A below) . its first partial 
derivative with respect to r' ■ is the radial slope 
(Equation 21B below). Its second partial derivative 
5 with respect to r' is the radial curvature 

(Equations 20C above) . 

spline' (r' ,6) = Co (9) + ci(0)r f + c 2 (9)r' 2 + C3 (9)r' 3 

Equation 21A 

d spline' (r' , 9) /dr' = C i(9)+ 2c 2 (9)r' + 3c 3 (9)r' 2 
0 Equation 21B 

The spline coefficients c n (9) are solved at each 
meridian 9 with the boundary conditions at the edge of 
the ablation zone, r' = R AZ ' and at the edge of the 
optical zone r' = R oz ' 
5 c 0 (9) + C!(9)R AZ ' + c 2 (9)R AZ ' 2 + c 3 (9)R AZ ' 3 = 0 

Equation 22A 
Ci(6)+ 2 c 2 (9)R A2 ' + 3 c 3 (9)R A2 ' 2 = 0 

Equation 22B 

c 0 (9) 4- Ci(9)Roz f + c 2 (6)R 0 z /2 + c 3 (9)R 0 z' 3 = a' (Roz' , 6) 
0 Equation 23A 

ci(9)+ 2 c 2 (9)Roz' + 3 c 3 (9)R 0 z' 2 = b' (R oz r , 9) 

Equation 23B 

where 
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b' (r',0) = 5a' (r',9)/dr' is the radial slope of the 
ablation profile in local polar coordinates. 

Equations 22 and 23 form a system of 4 linear 
algebraic equations with 4 unknowns c n (e). These 
5 equations are solved using well known linear algebra 

techniques . 

Depth Constant; :for Non- Posit:! vity Constraint 

The ablation depth, in accordance with the present 
convention, is non-positive. A negative number 

10 . corresponds to tissue removal. A positive number 

corresponds to tissue addition, which is not possible 
with laser ablation. To enforce this non-positivity 
constraint, a constant depth offset d is added to the 
entire ablation map. d is computed in two steps. 

15 d = di + d 2 

Equation 24 

where 

di + d 2 are the two components of d. 

di = -maximum{ a 0 ' (r',9)} for r' < R 0 z' 
20 Equation 25A 



ai' (r',9) = a 0 ' (r',9) + d x for r' < R 02 ' 

Equation 25B 



where 
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maximum! } is a function that returns the maximum 
value of the enclosed function, and 

ao' (r',9) is the ablation map without any constant 
offset. 

5 Equation 25A computes the minimum depth offset to 

keep the OZ non-positive. 

d 2 = -maximum! (1/3) bi'(Roz',9) W TZ ' + ai'(R 0 z',9)} 

Equation 26A 
a' (r',9) = ai'(r',9) + d 2 for r' < R oz ' 

0 Equation 26B 

where 

bi'(r',9) = 5ai' (r' ,9) /dr' is the radial slope of 
the ablation profile in local polar coordinates. 

Equation 2 6A computes the minimum additional depth 
5 offset required to keep the TZ non-positive when there 

is an upward radial slope at the edge of the OZ . 

Iterative Deconvolution Loop 

Iterative deconvolution is used to generate the 
ablation map. The ablation map is initially set to be 
0 the same as the correction map inside the OZ. The 

correction map Ah target (r , 9) specifies the desired 
corneal surface change (Equation 13) . It is converted 
to the local-coordinate form for the deconvolution 
computations . 
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a 0 '[0](r',9) = Ah' target (r, 6) for r' < R 02 ' 

Equation 27 

where 

a 0 ' [0] (r',9) is the initial ablation map (the 
5 bracket [0] denotes the 0 th iteration) . 

The depth offset is then computed according to 
Equations 25A-26B and the transition zone map is 
computed according to Equations 19A-23B to yield the 
initial ablation map a' [0] <r' , 9) . 
10 The corneal surface smoothing is then simulated by 

convolution (Equation 14) . This is done in Cartesian 
coordinates . 

Ah']ilx\,y') = aii-llx\y)<2>f'(x\y) 

Equation 28 

15 where the iteration index bracket fi] denote the 

ith iteration. The loop starts with i = 1. 

The smoothing effect is canceled out by adjusting 
the ablation map. This is done in polar coordinates. 
g'[i](r',G) = Ah' [i] (r',9) - Ah' target (r ' , 9) for r' < R 02 ' 
20 Equation 29A 

a 0 ' [i + 1] <r' ,9) = a 0 ' [i] ( (r' , 9) - g' [i] (r' , 9) 

for r' < R 0 z' 

Equation 29B 
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where g' [i] ( r ' ,9) is the difference map between the 
outcome map Ah' [i] (r',9) and the target correction map 

Ah' target (r' ,9) . 

Again, the depth offset is then computed according 
5 to Equations 25A-26B and the transition zone map is 

computed according to Equations 19A-23B to yield the 
compute initial ablation map a' [i+1] (r',9). The loop 
is then repeated and the iteration index is 
incremented. The loop ends when the outcome and target 
corrections have the same shape and the difference map 
approaches a constant value. 

5 = maximumlg' [i] (r' , 9) } - minimum{g' [i] (r' ,9) } 

for r' < R oz ' 

Equation 30A 

15 5 < 6 max 

Inequality 30B 
Inequality 30B is the termination condition for 
the iteration loop. The maximum tolerable difference 
5max is a small depth, such as 0.1 jam, that is smaller 
than the resolution of the ablation process. The final 
iteration yields an ablation map a' (r'0) that is very 
close to the ideal needed to produce target correction. 
The ablation map is converted to the global Cartesian 



20 
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coordinate (Equations 15A-16D) for the next step, where 
the laser pulse sequence is calculated. 

The above algorithm is used to generate ablation 
maps for myopia (Fig. 4A) , hyperopia (Fig. 4B) , myopic 
5 astigmatism (Fig. 4C) and hyperopic astigmatism 

(Fig. 4D) . The smoothing constant s = 0.5 mm was used. 
For myopia, the deconvolution deepens the ablation in 
the OZ. This is caused by steepening of the ablation 
slope in the peripheral OZ (Fig 2A, 4A) . The 

10 algorithm, in accordance with the present invention, 

compensates for the healing response after any 
correction pattern, including myopia, hyperopia, 
astigmatism, and higher-order aberration. 

For hyperopia correction, the deconvolution also 

15 deepens the ablation in the OZ, more in the periphery 

than at the center (Fig 2B, 4B) . The hyperopic 
ablation map, in accordance with the present invention, 
puts the deepest ablation outside of the OZ . This is 
significantly different from conventional ablation 

20 algorithms, which place the deepest ablation at the 

edge of the OZ . The deconvolution algorithm, in 
accordance with the present invention, guarantees 
accurate full correction in the OZ and extends some 
correction into the transition zone. In contrast, 
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10 



conventional hyperopic ablation achieves no correction 
in the TZ and produces full correction only in the 
center of the OZ, with gradual decrease in correction 
toward the peripheral OZ . 

To achieve the improvement in optical quality in 
the OZ, the ablation patterns, in accordance with the 
present invention, require more ablation depth than 
conventional algorithms. The ablation depths per 
diopter for the correction of myopia and hyperopia are 
tabulated for a range of OZ diameters in Table 1. The 
depth is measured at the deepest point of the ablation. 



15 



OZ size 


5mm 


5 . 5mm 


6mm 


6 . 5mm 


Myopia 


12.4 


14 . 4 


16.6 


18 . 8 


Hyperopia 


14 .4 


16. 6 


19.1 


NA 



20 



— * f ^- ^j-w^u^i. , xucaouicu cl l tne 

deepest point of the ablation map. 

The myopic, astigmatism ablation map in accordance 
with the present invention has an elliptical AZ 
(Fig. 4C) and a circular OZ . The eccentricity of the 
ellipse is determined by the ratio of astigmatism to 
myopia (Equations 20) . The ablation profile on the 
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steep meridian is similar to that for pure myopia 
(Fig. 3B) . On the flat meridian, the ablation is 
deeper at the edge of the OZ to compensate for the 
expected epithelial hyperplasia at that location 
5 (Fig 3A) . 

The hyperopic astigmatism ablation profile has the 
same shape as that for hyperopia, but the depth of 
ablation varies from meridian to meridian (Fig. 4D) . 

Pulse Sequence Generation 

10 Laser systems produce ablation patterns by 

scanning laser pulses on the cornea. Thus, the final 
instruction to the laser system consists of a sequence 
of locations for the placement of pulses. In 
accordance with the present invention, the ablation 

15 maps are translated into pulse placement sequences with 

a "rain-drop" algorithm. An optional sprting algorithm 
is then used to order the pulses so that temporally 
nearby pulses do not overlap spatially. 

The simplest method of translating an ablation map 

20 to a pulse map is to apply pulses with a density 

proportional to the ablation map a(x,y). This approach 
leads to a distortion of the ablation caused by the 
smearing effect of the laser spot. The larger the 
spots used, the more the outcome deviates from the 
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target ablation. The distorted outcome can be 
simulated by convolving the target ablation map with 
the laser spot map. To reduce the distortion, a ^rain- 
drop" .algorithm is used to plan the pulse map. The 
5 rain-drop algorithm is a constrained iterative 

deconvolution procedure, like the algorithm used to 
generate the ablation map. A major difference is that 
the pulse map necessarily consists of discrete impulses 
to represent the pulses, whereas the ablation map is 
10 continuous. Another difference is that the pulse map 

is linked to a pulse sequence that has an additional 
temporal dimension. The rain-drop algorithm introduces 
a random element to give the sequence a desirable 
fractal property. 

15 Rain-Drop Algorithm 

The rain-drop algorithm generates randomly located 

pulses. Whether the pulses are rejected or lands on 

the simulation grid depends on several constraints and 

a probability test. The constraints are: 

20 1- Non-coincidence. Pulses cannot fall on 

the exact same location twice. This 
condition is not essential, but helps 
reduce deep edges, which is an issue 
with flat-top beam profiles. 

25 . 

2. Non-positivity . Pulses cannot fall 

where there is positive value within the 
spot area. 
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Border. The entire spot must fit inside 
the AZ borders. 



If a pulse passes the constraints, a probability test 
5 is applied to determine if the pulse. is finally 

accepted. The probability of acceptance is 
proportional to the average depth of the remainder map 
inside the pulse spot. The averaging may be weighed by 
the pulse map (more accurate) or weighed uniformly 

10 (faster computation) . The remainder map is the 

difference between the achieved .ablation and the target 
ablation map. 

The rain-drop algorithm can be applied to a 
situation where more than one spot size is used. The 

15 largest spot is used first. The sequence switches to 

the next smaller spot size after a preset number of 
larger spots are rejected by the constraints and 
probability test. After the random rain-drop algorithm 
runs through all the spot sizes, a deterministic fill- 

20 in phase is used to fill in the remaining low spots on 

the remainder map. The fill-in algorithm only uses the 
smallest spot size. The constraints are not applied in 
the fill-in phase. The remainder map is convolved with 
the spot map to obtain a guide map. Pulses are applied 
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to the lowest spot on the guide map until the average 
level of the remainder map reaches zero. 

The random nature of the rain-drop algorithm gives 
the generated pulse sequence a fractal distribution. 
5 Fractal patterns are self-similar at different size 

scales. Fractal distribution means that the pulse 
sequence is ordered in such a way so that if the 
sequence is stopped at any point, the fractional 
correction has approximately the same shape as the 
complete target correction. This property also makes 
the shape of the overall correction less susceptible to 
temporal variations in pulse energy, tissue hydration, 
or any other aspect of the ablation process. 

The rain-drop algorithm is detailed in the 
15 following program code written in the style of C 

programming language . 

r RAINDROP ALGORITHM TO GENERATE PULSE MAP & SEQUENCE V 

' nt i> ?* indices for sequence & loops */ 

int ix » j y; /* indices for map arrays 7 

2 0 int ix_start, iy_start; r starling indices for the ablation zone 7 

float target_volume. meanjeve!, bottorn_depth, az_area, az_width, az_height; 

float spot_diameter[N_SPOT_SIZES]; r spot sizes (diameters) 7 

float spoLvolume[N_SPOT_SIZES]; r ablation volume per pulse for each spot size 7 

float spot_mapfN_SPOT_S!ZESI[NX][NY3; r ablation maps for a single pulse 7 

2 5 /• NX, and NY, are the number of elements in the x and y dimensions of the simulation grid V 

float pu!se_map[N_SPOT_S!ZES] [NX][NY]; r one pulse map for each spot size V 

float remainder[NX][NY]; r remainder map V 

float guide_map[NX][NY]; 
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float target_ablation[NX][NYl; r the total target ablation map */ 

struct Location Jndices {int ix, iy;} r structure consisting of x, y indices V 
Locationjndices minjxy /• location indices of the minimum value on the map V 

struct Location {float x, y;} r structure consisting of x, y coordinates V 

5 Location sequence[N_SPOT_SIZES][MAX_PULSES] r Pulse sequences 7 

int N[N_SPOT_SIZES]; r the number of pulses for each fraction of pulse map */ 

int rejection, rejection_flag, max_rejection; 
r initially, remainder map = ablation map */ 
equate_map(remainder, target_ablation); 
10 r equate_map(map1 , map2) performs the array function mapl = map2 */ 

target_volume = Integra te_volume(target_ablation); 
r calculate the volume (negative) of the remainder map 7 

meanjevel = targe t_volume / az_area; /* az_area is the area of the ablation zone */ 
bottom_depth = min_map( remainder); r returns the most negative value on remainder map V 

15 

f* random "rain-drop" stage */ 

for (i=0; i < N_SPOT_SI2ES-1; i++) { 

make_spot_map(spot_mapp], spot_diameter[i]); r generate spot map V 

spot_volume[i] = integrate_volume(spot_mapti)); r spot volume is negative V 

2 0 j = 0; f* sequence count reset to zero V 

rejection = 0; /* consecutive rejection count V 

max_rejection = az_area/(Pr(spot_diameter[i]/2) A 2) + 100; 

while (rejection < max_rejection) { 
2 5 r randomO generate random value between 0 to 1 7 

r DX and DY are grid element sizes along horizontal and vertical dimensions */ 

r ix_start and iy_start are the lowest x.y indices where the A2 starts */ 

r az_width and az_height are the horizontal and vertical diameters of the AZ */ 

ix = ix_start + round(random()*az_width/DX); r round(float) returns the closest integer 

30 v 

'y = iy^start + round (random()*az_height/DY); 
rejection_flag = 1 ; 
if (pulse_map[i]fix][iy] == 0) 
r pulse cannot fall on same place twice */ 
35 if (inside_AZ(ix,iy,spot_diameter[i] 12) == TRUE) { 

I* pulse must not ablate outside the ablation zone V 
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if (spot_max(remainder, ix, iy, spot_diameter[i]) < 0) 

/* spot_max() return the maximum value inside the pulse spot diameter. 

Pulse is rejected if remainder level above zero anywhere inside spot. 7 
if (spot_mean(remainder,ix,iy,spot_diameter[i]) / bottom_depth > random()) { 
5 r spot_mean() calculate the mean depth inside the pulse spot diameter V 

sequence[i]0++] = locate(ix, iy); 

r the location x.y conesponding to indice ix.iy is added to the pulse sequence 7 

puise_map[i][ix][iy) = 1; 

ablate(remainder, spot_mapp], ix, iy); 
^ 0 r ablate(Map.SpotJx.iy) subtract Spot at Locationjndices ix.iy from Map 7 

mean Jevel -= spot_volume[i] / az_area; 

bottom_depth = mi n_map( remainder); 

rejection_flag = 0; 

rejection = 0; 
15 } 

if (rejectionjtag = 1) rejection ++; 

' » 

N[i]=j; 
} 



20 



/• deterministic fill-in stage for the last layer of the last spot size (smallest) 7 



while (meanjevel < 0) { 

guide_map = convolve(spot_map[i], remainder); 

2 5 minjxy = minJocation_indices(guide_map); 

r minjocationjndices(map) returns the x,y indices of the minimum value on the map 7 
sequence[i][j++] = locate(minJxy.ix f minjxy.iy); 
pulse_map[i][min_ixy.ix]fmin_ixy.iy] += 1; 
ablate(remainder, spot_map[i] t min_xy.ix, min_xy.iy); 

3 0 meanjevel -= spot_volume[i] / az_area; 

} 

ND] = j; 

r END ALGORITHM V 

Because the rain-drop algorithm constantly updates the 
35 remainder map to guide subsequent pulse placement, the 
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final remainder map is very close to flat zero. The 
accuracy of the algorithm is demonstrated using the 
following examples . 

For the examples, a Gaussian beam v fluence profile 
5 is used defined in the following equation. 

G(r) = F 0 exp{-r 2 /Rbeain 2 } 

Equation 31 

where 

^beam is the 1/e beam radius, 
10 exp( ) is the natural exponential function 

with base e, and 

F 0 is the peak fluence. 
The ablation profile is calculated from the beam 
fluence using a Beer's law approximation of the 
15 ablation process. The ablations characteristics of 

the 193 nm wavelength ArF excimer laser were used. 

S(x,y) = E ln{G(r)/F th ) ) for G(r) > F th 

Equation 32 
S (x, y) = 0 for G (r) < F th 

20 where 

S (x,y) is the ablation depth map from a single 
laser pulse ("spot map"), 

ln{} is the natural logarithm function with 
base e, 
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E is the ablation efficiency, and 

F th is the threshold fluence for ablation. 

The spot diameter is defined as the diameter 
inside which S{x,y) > 0. 
5 For our simulations, an ablation efficiency of 0.3 

micron and a threshold fluence of 60 mJ/cm 2 are adopted. 

For the examples, a peak fluence F 0 = e F th was 
used. This provides a central ablation depth of 0.3 
microns and ablation spot radius R spot = R beam . Ablation 
0 spot diameters of 2 . 0 mm and 1 . 0 mm were used. The 

rationale for using two spot diameters is to use the 
larger spot to remove most of the required volume 
quickly, then use the smaller spot to achieve a higher 
final accuracy. 

5 The pulse maps were generated for 1.0 D of myopia 

(Fig. 5A) and hyperopia. The rain-drop algorithm was 
able to use the larger 2-mm spot for more than half of 
the pulse counts and for' the great majority of ablation 
volume (Table 2) . 
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Spot 

Diameter 


2 . 0 mm 


1 . 0 mm 


Total 


myopia 
(-1.00 D) 


# spots 


524 


4 56 


980 


% volume 


82. 1 


17 . 9 


100 


hyperopia 
( + 1.00 D) 


# spots 


872 


850 


1722 


% volume 


80.4 


19.6 


100 



Table 2. The distribution of pulses generated by the rain- 
drop algorithm using two different spots diameters. See Figs. 4, 
5 for corresponding ablation and pulse maps. 



5 

The remainder maps for the ID ablations show 
variations on the scale of the 0.3 micron pulse 
ablation depth. The myopic remainder map is shown in 
Fig. 5B. The remainder map is the difference between 

10 the achieved ablation (pulse map convolved with the 

spot map) and the target ablation map a(x,y). The root- 
mean-square {RMS ) values of the remainder maps for both 
the myopic and hyperopic examples are below the 0.3 
micron pulse ablation depth (Table 3) . After the 

15 expected surface smoothing, the achieved corrections 

deviate very little from the target corrections. The 
smoothed residual map for the myopic correction (Fig. 
5C) shows a very slight under-correction pattern. The 
under-correction of 0.06 micron RMS is equivalent to - 

20 0.02D , which is clearly insignificant. The RMS 

smoothed residual for the hyperopic correction is 
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similarly very small (Table 3) . The smoothed residual 
map is defined as the difference between the achieved 
correction (pulse map convolved with the spot map and 
then smoothed according to Equation 14) and the target 
5 correction map (Equation 13) . 



myopia 
(-1.00 D) 


Raw 


0.136 


Smoothed 


0.058 


hyperopia 
(+1.00 D) 


Raw 


0.131 


Smoothed 


0.062 



Table 3. Root-mean-square (RMS) residual in simulated 
ablations using 2 & 1 nun diameter spots. 



Sequential Non-Overlap Sorting Algorithm 

10 It is also desirable to order the pulse sequences 

so that consecutive pulses are placed in non- 
overlapping locations. This prevents the plume of a 
laser pulse from affecting the transmission of 
subsequent pulses. The rain-drop algorithm already 

15 yields sequences with little sequential overlap. 

However, an additional sorting algorithm is preferably 
used to further minimize sequential overlap. This 
sorting operation is detailed in the following program 
code. First, the number of pulse intervals where the 

20 non-overlap condition must be observed is calculated by 

dividing the plume dispersal time by the laser firing 
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interval and rounding up the result to the next larger 
integer. The pulses are examined consecutively in 
reverse order. The order was reversed because the 
clumping of pulses tends to occur late in the rain-drop 
5 algorithm. The pulse being examined, is the index 

pulse. The spatial distance between the index pulse 
and the next pulse is compared with the required 
minimum separation (the spot diameter or slightly 
larger) . If the separation is too small, the next 
10 - pulse is exchanged with the nearest later pulse that 
satisfied the required separation with the index and 
previous pulses . The index is incremented until the 
end of the sequence is reached. 

r ALGORITHM TO SORT PULSE SEQUENCE FOR CONSECUTIVE NON-OVERLAP •/ 
15 int i, j, k, m, t, flag, n_interval; 

Location temp; 

for (i = N_SPOT_S!ZES-1; i >=0; i-) { 

njnterval = ceiling(plume_dispersal_time[ij / PULSEJNTERVAL); 
min_separation = PLUME_SPREAD * spot_diameter[i]; 
20 forG = ND>1;j>0;H 

for (k = 0; k < njnterval; k++) 

if (distance(sequencep]GL sequence[i]G-k-1]) < min_separation) 
for (m = j-k-2; m >= 0, m++){ 

flag = 1 ; r If flag = 1 , swappable*/ 
25 for (t = j-k, min(N(i), j-k-1 +n_interval), t++) 

if (distance(sequence [i][t] r sequence 
[i][m]) <= min_separation) 
flag = 0; 

ifflag = 1{ 
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temp = sequencep]D-k-1]; 
sequencep][j-k-1J = sequence[i][ m]; 
sequencep][m] = temp; 
break; 

5 } 

} 

} 

r END ALGORITHM 7 
r COMMENTS 

1 0 plume_dispersa1_timep] is the time needed to disperse the plume so it does not interfere with later pulses, 

PLUME_SPREAD is the ratio between the effective plume diameter and the ablation spot diameter, 
njnterval is the number of pulse interval that must elapse before overlap is allowed, 
ceiling(number) returns the number rounded to the next larger integer, and 
distance(pulse1 , puise2) returns the distance between pulsel and puise2. 

15 v 

As an example, the sequence of the 2 mm spots in 
the 1 D myopia example are sorted. The non-overlap 
interval number was set to 1 (n__interval = 1) . The 
plume size was assumed to be the same as the ablation 

20 spot size ( PLUME_SPREAD = 1) . The mean distance between 

consecutive pulses before sorting was 2.70mm. Sorting 
improved this to 3.02mm. The mean area overlap between 
consecutive pulses was 7.25% before sorting. This 
improved to 0.27% after sorting. The sorting algorithm 

25 is highly effective in reducing overlap. It has 

minimal effect on the fractal distribution of the pulse 
sequence because the swaps are few and random. 

Referring to Figure 6, a system is shown for 
practicing the present invention. A laser source 10 
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genejates laser pulses that travels along beam path 12. 
These laser pulses are of the appropriate wavelength, 
duration, and energy for the intended corneal ablation. 
For example, an argon fluoride excimer laser operating 
5 at 193 nm can be used, with pulse duration in the 

millisecond range. The beam path is directed by mirror 
or mirrors 20 to beam shaping apparatus 30. Beam 
shaping is performed by a combination of an aperture 
32, and lenses 34 and 36. The beam is then .steered by 

10 scanning apparatus 40. The scanning apparatus 

preferably consists of two steering mirrors 42 and 44 
along two perpendicular axes. The beam is then 
•deflected by wavelength-selective mirror 50 onto the 
target eye 60. During the operation of the laser, the 

15 position of the eye is monitored by aiming apparatus 

70. The aiming apparatus preferably consists of a 
microscope through which the surgeon can visualize the 
target eye, a video camera and other tracking sensors. 
The position of the eye can be measured by computer 

20 processing of the video images, or other types of 

tracking sensor output. Computer 80 is operatively 
connected to the beam shaping apparatus 30 and the 
scanning apparatus 40 for controlling those devices. 
The output of the arming apparatus provides sensor 
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information to computer 80. Computer 80 is operated by 
the laser operator through an appropriate input 
device 82 under the direction of the surgeon. The 
computer stores a pulse sequence and control algorithms 
5 described above which are used to control the laser- 

pulse energy and timing, beam shaping, and beam 
steering. Input from aiming apparatus 70 is used to 
yield tracking information that modifies the 
instruction sent by computer 80 to steering 

10 apparatus 40. 

An eye 60 is illustrated in Figure 7. The 
cornea 90 is the target of the laser surgery used to 
correct the refractive error of the eye. The line of 
sight is the line connecting the center of the pupil 

15 with the object of regard. The pupil is the aperture 

formed by the iris 92. The Radius r of the global 
coordinate system is defined by the perpendicular 
distance (radius of the pupil) from the line of sight. 
The radius r' of the local coordinate system is defined 

20 as the arc length from the line of sight along the 

anterior corneal surface. The z-axis of the global 
coordinate system is parallel to the line of sight. 
The z' -axis of the local cornice system is defined as 
being perpendicular to the corneal surface. 
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Figure 8 is a top level flow chart of the ablation 
process 100 in accordance with the present invention. 
Conventionally clinical manifest refraction 102 is used 
to determine the setting of laser vision correction. 
5 Manifest refraction is determined by placing trial 

lenses in front of the. subject eye and changing the 
lenses based on the subject's evaluation of the 
resulting vision, until the best combination of lenses 

r 

is found. This combination of spherical and 
10 cylindrical lenses is termed the manifest refraction 

prescription. Manifest refraction can be performed 
after the eye ciliary muscle has been paralyzed with 
eye drops. This is termed "cycloplegic refraction". 
The surgeon may also modify the laser setting using 
15 keratometry 104, which measures the corneal curvature 

at one single diameter. Wavefront sensors 106 can be 
used to measure the detailed optical aberration of the 
eye, including spherocylindrical refraction and higher 
order aberration. The output of the wavefront 
20 sensor 106 is a map of wavefront height, which can be 

further analyzed to provide the amplitude coefficients 
for Zernike series. Corneal topography 108 provides a 
map of the anterior corneal surface height. This 
provides more detailed information than keratometry and 



10 



15 



20 
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can also be analyzed to provide Zernike coefficients. 
All this measured and derived data can be combined to 
provide a prescription for laser correction. If the 
prescription only consists of spherocylindrical terms 
(or the equivelant Zernike 2 nd order ;. terms ) then it is a 
traditional refractive prescription. If it contains 
higher order Zernike terms, this is a wavefront 
prescription. All prescription information is input to 
a nomogram algorithm 110 which is a correcting 
algorithm used to calculate actual settings for a laser 
system to achieve correction. The output 112 of the 
nomogram algorithm is referred to as the correction 
prescription. The correction prescription 112 is the 
input for the correction map-generating algorithm 114 
of the present invention. The resulting output 116 is 
the target correction map, which is specified in terms 
of the desired corneal surface height change. The 
target correction map output 116 is input to an 
ablation map-generating algorithm 118 that converts 
this to a target ablation map 120. The target ablation 
map 120 specifies the depth of laser. ablation to be 
performed. 

Referring to Fig. 9, details of the ablation map 
generating step 118 are shown in detail. The target 
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correction map output 116 is provided to a conversion 
step 130 that converts the map to local coordinates and 
provides the -target correction map in local coordinates 
to a loop initialization step 132. In step 132, the 
5 ablation map value is set equal to the target 

correction map and the difference map value is set 
equal to zero. The ablation map 133 is the stored 
ablation map that is being operated on by the' 
subroutine at the various steps. Those skilled in the 
.10 art will appreciate that the ablation map is called up 

from storage, operated on in accordance with a process 
step and then returned to the same storage memory. In 
step 134, optical zone construction occurs in which the 
ablation map value is set equal to the ablation map 

15 value minus the difference map value. This provides an 

ablation map without a transition zone. In step 136, 
the transition zone is constructed using cubic spline 
curves discussed above which provide smoothing. Next, 
a convolve ablation map with epithelial smoothing 

20 function is performed in step 138 to derive an achieved 

correction map. The difference map value is then set 
equal to the achieved correction map minus the target 
correction map in step 140 that is a difference map 
with a transition zone. The region outside of the 
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optical zone is cropped in step 142 to yield a 
difference map without a transition zone. In step 144, 
a determination is made as to whether the difference 
map range is less than a predetermined tolerance value. 
5 If the determination is negative, the process loops 

back to step 134. If the determination is affirmative, 
the process then converts the map back to global 
coordinates in step 14 6 and outputs the target ablation 
map 120. 

Referring back to Fig. 8, the target ablation 
map 120 is input to a pulse sequence generating 
step 150 that converts the target ablation map into a 
pulse sequence 152. The pulse sequence 152 specifies 
the location, diameter, and fluence profile of laser 
15 pulses. In addition, the pulse sequence contains 

order, size, and locations of the laser pulses. 

Referring to Fig. 10, details of the pulse 
sequence generation step 150 is shown. As discussed 
above, the largest spot size of the target ablation map 
120 is first analyzed using the rain-drop algorithm 
160. In the rain-drop algorithm, the rejection count 
value is set equal to zero in step 162. Randomly 
placed pulses are then generated in step 164. A 
determination is made in step 166 as to whether the 
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pulse meets constrains. The constraints are (1) 
Ablation Zone--the pulse fits within the ablation zone, 
(2) Non-Coincidence — the pulse does not coincide with 
other pulses, and ( 3 ) Non-Positivity — the - pulse does not 
5 cause any point on the remainder map to become 

positive. Recall, only negative values are permitted. 
If the determination in step 166 is .affirmative, a 
determination is made in step 168 as to whether the 
probability of acceptance is greater than random. The 

10 probability of acceptance is the average depth of the 

remainder map inside the pulse spot divided by the 
deepest depth on the entire remainder map. The random 
number is chosen from a uniform distribution on the 
interval between zero and one. 

15 If the determination is step 168 is affirmative, 

the process sets the rejection count to zero and adds 
the pulse value to the sequence in step 170 and the 
remainder map is updated in step 172. The remainder 
map is the difference between the achieved ablation and 

20 the target ablation map. If either determination 166 

or 168 is negative, the process increases the rejection 
count by 1 in step 174. From either step 172 or 
step 174, a determination is made in step 176 as to 
whether the rejection count has reached a predetermine 
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tolerance. If the determination in step 176 .is 
negative, the process loops back to step 164. If the 
determination in step 176 is affirmative, a 
determination is made in step 180 as to whether the 

5 smallest spot has been processed. If the determination 

is negative, the next smallest spot is determined in 
step 182 and the rain-drop algorithm is run on the next 
smallest spot. 

If the determination in step 182 is affirmative, 

0 the fill-in algorithm is run in step 190. In the fill- 

in algorithm 190, a convolved remainder map with spot 
map is provided in step 192. The remainder map is then 
updated in step 194 by placing the pulse at the minimum 
point on the index map. A pulse is then added to the 

5 sequence in step 196. A determination is then made in 

step 198 as to whether the average depth of the 
remainder map is positive. If the determination is 
negative, the process loops back to step 192. If the 
determination in step 198 is affirmative, the process 

0 then retrieves the stored pulse sequence in step 200 

and the pulse sequence 152 is then provided. 

Referring back to Fig. 8, the pulse sequence 152 
is provided to the sequential non-overlap sorting 
step 210. Referring to Fig. 11, the details of the 
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sequential non-overlap sorting step 210 will be 
appreciated. The order of the pulse sequence is 
reversed in step 220 and the smallest spot is 
identified in step 222 for further processing. In 
5 step 224, the non-overlap sorting algorithm is 

performed. In this sorting algorithm, the first pulse 
is used set as the index pulse in step 230. In step 
232, a determination is made as to whether there is an 
acceptable separation between the index pulse and the 

10 next pulse. If the determination in step 232 is 

negative, the process proceeds to step 234 where a 
later pulse with adequate separation from the index 
pulse is identified. A determination is then made in 
step 236 as to whether such a later pulse, i.e., the 

15 pulse searched for in step 234, has been found. If the 

determination in step 236 is affirmative, the process 
then swaps the qualified pulse with the pulse next to 
the index pulse in step 238. From either an 
affirmative determination is step 232, a negative 

20 determination in step 236, or from step 238, the 

process makes a determination in step 240 as to whether 
the end of the sequence for the current pulse size has 
been reached. If the determination is step 240 is 
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negative, the index pulse is advanced in step 242 and 
the process then loops back to step 232. 

If the determination in step 240 is affirmative, 
the process makes a determination in step 250 as to 
5 whether the spot size processed was the largest spot 

size. If the determination in step 250 is negative, 
the process goes to the next large spot size in 
step 252 and then . loops back to step 224. If the 
determination in step 250 is affirmative, the process 
10 then reverses the order of the pulse sequence in 

step 256 and outputs the final pulse sequence. 
Referring back to Fig. 8, the final pulse sequence is 
then provided to the laser control for vision 
correction . 

15 Construction of Mathematical Model to Account for 

Corneal Surface Smoothing After Laser Refractive Surgery 

A model of epithelial transport is proposed based 
on first principles. The first postulate is that 
20 epithelium migrates toward a depression on the corneal 

surface. Mathematically, this is similar to the flow 
of solute across a concentration gradient and can be 
described by the partial differential equation below. 
For a corneal surface region S bound by a closed line 
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L, the rate at which epithelial volume flows out of a 
line element dl depends on the surface height gradient. 

-dQ m = -^i (Vh' ) * n dl < 

where 

5 Q m is the migratory epithelial volume flow, 

ji is the motility of epithelium, 
V is the gradient vector operator in two 

dimensions ± d/dx' + j d/dy' , where .x' and y' are local 

transverse coordinates tangential . to the surface, 
10 h' is the corneal surface height along the z' 

axis, which is defined locally as normal to the 

surface, and 

n is the outward vector normal to the line 

element dl . 

15 The net migratory flow of epithelial volume into 

the area enclosed by boundary L is given by a line 
integral . 

Qm = J L V ( v ^' ) • n dl 

Equation 33 

Application of the divergence theorem in two 
dimensions gives us a surface integral. 

Qm = JJs M- (V 2 h' ) da 

Equation 34 
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where 

V 2 is the Laplacian operator in two dimensions 
cf/dx' 2 + d 2 /dy' 2 , and da is the surface element. 

Epithelial cells divide and provide replacement t 
5 the epithelial volume. It is further postulated that 

cell division occurs at a constant rate. Although thi 
may not be true in an active wound healing phase, it 
should be close to reality in the state of equilibrium 
In integral form, this event can be expressed as: 
10 Q g = fJs Y da 

Equation 35 

where 

Q g is the generative epithelial volume flow, and 
y is the generativity of epithelium (rate of 
15 epithelial growth or replacement) . 

Epithelial cells mature and move to the surface, 
where they eventually slough off. It is postulated 
that the rate of loss is proportional to the thickness 
of the epithelium p' . 

20 Qi = Us -A. p' da 

Equation 36 

where 

Qi is the epithelial volume flow due to loss, 



WO 03/075778 PCT/US03/06343 

-67- 

X is the lossivity of epithelium (rate of 
epithelial sloughing) , and 

p' ,is the epithelial thickness measured along the 
surface normal axis z' 
5 To complete the picture of epithelial transport, 

the epithelial thickness is linked to corneal surface 
height . 

h' = s' + p' 

Equation 37 

10 Q e = JJ S dp' /dt da • 

Equation 38 

where 

s' is the subepithelial surface elevation, and 
Q e is the absorption of epithelial volume flow by 
15 epithelial thickness change. 

Conservation of epithelium provides the master 
equation . 

Qe = Qm + Qg + Ql 

Equation 39 

20 -In integral form, the equation expands to: 

fj s dp' /dt - |i (V L 2 h') - y+ X p' dc = 0 

Equation 40 
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Since the model applies to any region S on the corneal 
surface, the equation can be written in the 
differential form where the integrand vanishes. 

dp' /St - n (V 2 h' ) - y + X p' = 0 . 
5 Equation 41 

In equilibrium, dp' /dt = 0 and therefore 
- \i (V 2 h' ) - y+ k p' = 0 

Equation 42 

The equation is simplified to obtain measurable 
10 constants. 

p' = c + s 2 V 2 h' , 

Equation 43 

where 

c = y/X represents the constant component of the 
15 epithelial thickness determined by a balance between 

growth and loss. The constant c has a unit of length 
(thickness) . 

s = V(|i/A.) is the smoothing constant. It has a 
unit of length and can be thought of as the radius over 
20 which smoothing occurs. It is .determined by a balance 

between epithelial migration and loss. One can also 
think of it as the typical distance over which 
epithelium migraites before it sloughs off. The term s 2 
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V 2 h' describes the variable component of the epithelial 
thickness that responds to surface curvature changes. 

The above model is applied to the prediction of 
corneal surface smoothing after LAS IK and PRK. 
5 Equation 44 below states that the post-operative 

corneal surface height hi(x',y') is equal to the pre- 
operative height h' 0 (x' , y' ) plus the ablation a' (x',y' ) 
(depth as negative quantity) and the change in 
epithelial thickness . Ap' (x' , y' ) . 
10 h'i = h' 0 + a' + Ap' 

Equation 44 

where Ap' = Pi' - po' 

Using Equation 43, we find that the change in 
epithelial thickness is: 
15 Ap' = s 2 V 2 Ah' 

Equation 45 

where 

Ah' =hi' - h 0 ' 

Applying Equation 4 5 to Equations 44, we find 
20 Ah' = a' + s 2 V 2 Ah' 

Equation 4 6 
Equation 46 has solutions in the forms of 
separable complex exponential functions exp(jco y /x') 
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= cos(co x 'x') + j sin(co x 'x) and exp(jco/y') = cos (go/ y' ) 
+ j sin (©/ y'). Therefore, it is useful to look at it 
in the frequency domain. Equation 4 7 is derived from 
Equation 4 6 using the differential property of the 
5 Fourier Transform. 

AH' = A' + s 2 (CD/ 2 + (D/ 2 ) AH' 

Equation 4 7 

where 

AH' (co x ' , cd/ ) is the 2-dimensional Fourier 
10 Transform of Ah'(x',y'), 

A' (cd/ , cd/.) is the 2-dimensional Fourier 
Transform of a'(x',y'), and 

co/ and cd/ are the respective spatial radian 
frequencies for x' and y' in radians/length. 
15 Equation 47 can be expressed using a transfer 

function F' ((D x ' r co/ ) . s 

AH' = F' A' 

Equation 48A 

F' <©x' , 0) y ') = 1 / [1+ (0)//C0 c ' ) 2 + <CD/ /G) C ') 2 ] 

20 Equation 48B 

where the cutoff radian frequency is co c ' = 1/s in 
radian/length . 
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The transfer function F' is identical to the 
square of a first-order Butterworth low-pass filter in 
two dimensions. Thus, the action of the epithelium is 
theoretically identical to a low-pass filter. In other 
5 words, the low-pass filter F' characterizes the 

frequency response of the corneal surface. 

The local coordinate system x' , y' , z' has been 
used. The local coordinate system has the perspective 
of a small . epithelial cell on the corneal surface and 

10 is important to the development of the epithelial 

smoothing theory. However, the corneal height and 
ablation depth are usually specified in fixed global 
coordinate system x, y, z. In the global coordinate 
system, z is defined as the line of sight and x, y are 

15 perpendicular to z. The global coordinate system has 

the perspective of the laser scanning system looking 
down at the cornea. The coordinate systems are 
identical at the corneal apex but the deviation 
increases as the radial distance approaches the corneal 

20 radius of curvature. 

The following approximate transform equations are 
used to convert between the local and global 
coordinates. The transform treats the corneal surface 
as a perturbation from a sphere. These are 
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approximations given that the epithelial smoothing 
change is a small correction on the effect of the 
ablation and the transform is a secondary correction on 
the epithelial smoothing calculation. 
5 a'(r',0) = a(r,9) cos [arc'sin (r /R c ) ] 

Equation 4 9A 

r' = R c arcsin(r/R c ) 

Equation 49B 
Ah(r,0) - AhMr',6) / cos(r'/8c) 

Equation 49C 
r = R c sin (r' /R c ) 

Equation 4 9D 

where 

r is the radius perpendicular to the central axis 
15 defined by the line of sight, 

R c is the radius of curvature of the corneal best- 
fit sphere (a value of 7.6 mm is used in our 
simulations) , 

height h is measured along the z axis, which is 
positive outward along the line of sight, 

x is perpendicular to the line of sight and 
positive to the right of a viewer facing the eye, 
y is perpendicular to the line of sight and 
positive upward, 
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G is the meridian angle, which is zero on the +x 
axis and increases counterclockwise, and 

r' is the arc length from the corneal apex. 
The coordinate transforms in Equation 49 are 
5 performed in polar coordinates and the low-pass filter 

in Equation 48 is specified in Cartesian coordinates. 
To convert between the polar and Cartesian coordinate 
systems, the following equations are used. 

x = r cos 9 

10 Equation 50A 

y = r sin 0 



15 9 = arctan(y/x) 



x' = r' cos 9 



20 



y' = r' sin 9 



Equation 50B 



Equation 50 C 



Equation 50D 



Equation 51A 



Equation 51B 



Equation 51C 



WO 03/075778 PCT/US03/06343 

-74- 

9 = arctan (y' /x' ) 

Equation 51D 

Ablation simulations can be performed using MatLab 
software (The Mathworks, Inc. Natick, MA, USA) . The 
5 corneal surface is simulated as a sphere of 7.6 mm 

radius. Surfaces, ablations, and surface changes are 
computed on digital grids of 10 x 10 mm with a sampling 
interval of 0.02 mm. The ablation maps are based on 
modifications (see specifics below) of the exact 

10 Munnerlyn's algorithm for target 1 diopter (D) 

corrections. The ablation map is Fourier transformed 
to the frequency domain. The smoothing action is 
simulated with the low-pass filter in Equation 4 8 using 
smoothing constants ranging between 0 and 1 mm. Then 

15 the corneal surface height change in converted from the 

frequency domain to the spatial domain. The smoothed 
surface change is applied to the simulated corneal 
surface . 

The surgically-induced refractive change is 
20 evaluated using Zernike polynomial decomposition 

Zernike polynomial decomposition is applied to the 
corneal surface height to obtain coefficients h n , ra , 
which are the amplitudes of the Zernike terms z n m . The 
change coefficients &h n , a are computed by taking the 
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difference between post-operative and pre-operat ive 
coefficients. The achieved refractive effect of the 
corneal surface change is evaluated with the formulae 
given by Schwiegerling et al. (J Opt. Soc. Am. A, 1995, 
5 "Representation of Videokeratosdocopic Height Data With 

Zernike Polynomials'' 12(10):p. 2105-2113) with the 
substitution of corneal refractive index in place of 
keratometric index. The change in the coefficient Ah 2 ,o 
of the def ocus term Z 2 ° is used to evaluate spherical 

10 equivalent ( SE) change (Equation 52A) . The change in 

the coefficients Ah 2/±2 of .the astigmatism terms Z 2 ±2 is 
used to evaluate astigmatism correction (Equation 52B) . 
A 5.0 mm pupil diameter is used for the Zernike 
analysis of refraction because it is the average pupil 

15 size in the dim lighting from a projected eye chart. 

Five mm is also in between the average pupil sizes in 
darkness (5.8-6.1 mm) and in room light (4.1 mm). The 
correction/ablation ratio is the achieved dioptric 
correction for the ID target ablation. The ratio is 

20 normalized to 1 for s = 0. 



4J5 



Equation 52A 
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bAstigmatism = (n - 1) ^ ^(Ah 22 ) 2 _ 2 ) 2 

Equation 52B 

where 

ASE = change in spherical equivalent in dioptric 

5 unit, 

Aastigmatism = magnitude of astigmatism change 
vector in dioptric unit, 

R P = pupil radius in mm, a value of 5.0 is used in 
our calculations, 
0 n = corneal refractive index , a value of 1.377 

is used, 

Ah 2 , 0 = change in coefficient for Z 2 ° in |im, and 
Ah 2 ,± 2 = change in coefficient for Z 2 ±2 in |im. 
Spherical aberration is assessed. with the 
5 coefficient Ah 4 , 0 of the Zernike term Z 4 °. This was 

evaluated using both 5 and 6 mm diameter analytic zones 
(pupils) . 

Simulation of Hyper opic and 

Myopic Ablations With a IADARVislon 

0 Simulations of myopic (-1D) and hyperopic (+1D) 

ablations were performed for comparison to the observed 

degrees of under-corrections with the LADARVision laser 

system (Alcon Summit Autonomous, Orlando, FL) . The 
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hyperopic ablation has a 6.0 mm diameter optical 
zone (OZ) and 9.0 mm transition zone (TZ) . The exact 
Munnerlyn algorithm is used in the OZ . A cubic spline 
transition starts and ends with slope of zero in the T 
5 (Fig. 12) . The myopic ablation has a 6 . 0 mm diameter 

with no TZ (Fig. 13) . 

Simulation of Minus Cylinder Ablation Wirbh a EC-5000 

Simulation of 1 D of minus cylinder ablation on 
the EC-5000 laser (Nidek, Japan) was performed for 

10 comparison to published data. The EC-5000 ablation 

algorithm utilizes a circular OZ of 5 . 5 mm with a 
transition zone (TZ) out to 7 . 0 mm diameter. The lase 
ablation is shaped by both an expanding-slit and an 
expanding-circular aperture. The action of the slit 

15 aperture was simulated by a full Munnerlyn-algorithm 

cylindrical correction in the OZ with continuous 
tapering in the TZ . The circular aperture provides a 
constant-slope transition in the TZ . The ablation map 
is the product of the two aperture effects. The 

20 ablat ion profiles in the flat and steep meridians are 

.shown in Figures 14A and 14B, respectively. 

Clinical data analysis was based on a computer 
database for laser in situ keratomileusis (LASIK) . 
Cases were' reviewed over, approximately a 18 month 
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period for inclusion in the regression analysis. The 
inclusion criteria were: 

1. Spherical hyperopic or myopic 
corrections . 

2. No previous eye surgery. 

3. No astigmatism correction with the lase 
(refractive astigmatism < 0.5 D) . . 

\4. Pre- and post-operative best corrected 
visual acuity of 20/20 or better. 

5. Optical zone of 6.0 mm on laser setting 

6. No flap complication. 

The laser settings were based on manifest and 
cycloplegic refractions and corneal topography (C-Scan 
5 Technomed, Germany) . The Hansatome (Bausch & Lomb 

Surgical, St. Louis, MO) was used to create superiorly 
hinged flaps. Laser ablation was performed using the 
LADARVision laser system (Alcon Summit Autonomous, 
Inc., Orlando, FL) . Refraction obtained at the three 

10 month post-operative visit was used for outcome 

analysis. The database was maintained using spread 
sheet software. Statistical analysis was performed 
using JMP software Version 4 (SAS Institute, Cary, NC) 
.Surgically-induce refractive change (SIRC) was 

15 calculated by taking the difference between 

post-operative and pre-operative manifest refractions. 
The spherical equivalent (SE) component of refraction 
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was calculated by adding the spherical component to 
half of the cylinder magnitude. 

Results 

Refractive outcome for 19. cases of spherical 

\5 hyperopic LAS IK and 77 cases, of spherical myopic LASIK 

were analyzed. All cases utilized a 6 mm diameter OZ. 
Linear regression showed that each diopter of hyperopi 
ablation produced 0.708 D of effect (Fig. 15A) . The 
laser setting analyzed was the actual ablation setting 

10 that included laser's internal nomogram adjustment, 

which adds 50% ablation to entered hyperopic settings' 
up to 2 D and adds 1 D for entries of more than 2 D. 
Using the same laser, each diopter of myopic ablation 
produced 0.968 D of correction (Fig. 15B) . The 

15 clinical correction/ablation ratios of hyperopic and 

myopic ablations are significantly different (p. < 
0.001). Since the series are produced at the same 
center, with the same surgeon, laser, optical zone 
diameter, and siurgical technique, it is clear that 

20 hyperopic ablation is subject to more regression or 

other modifying effects than myopic ablation. 

' The simulations showed that the 
correction/ablation ratios for both hyperopic and 
myopic ablations decrease - with increasing smoothing 
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action, which is controlled in the model , in accordance 
with the present invention, by the smoothing constant 
s. With increasing s, the correction/ablation ratio 
fell faster for hyperopia (Fig, 16A) than for myopia 
5 (Fig. 16B) . But the difference was not. as great as the 

clinical data indicated (Fig. 4). A 
correction/ablation ratio of 0.708 for hyperopia 
correction is matched at s = 0.63 mm. A 
correction/ablation ratio of 0.968 for myopia 

10 correction is matched at s = 0.32 mm. 

The reduced correction is correlated with 
epithelial thickness modulation. In Figs. 12 and 13 , 
the epithelial thickness changes can be visualized as 
the difference between the achieved corneal surface 

15 change profiles and the ablation profiles. After 

hyperopic ablation (Fig. 12) f the model predicts 
epithelium thinning at the center and thickening at the 
periphery of the OZ, where the ablation is deepest. 
After myopic ablation (Fig. 13), the model predicts 

20 epithelium thickening at the center and thinning at the 

periphery of the OZ . 

The simulations show that, in the absence of 
smoothing (s = 0), there is induction of oblate 
spherical aberration (SA) after hyperopic ablation 
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(Fig. 17A) and prolate SA after myopic ablation 
(Fig. 17B) . This is expected from the spherical nature 
of Munnerlyn ablation algorithms. However, with 
increased smoothing, the opposite SA is induced. At 

.5 s = 0.5 mm, prolate SA is induced after hyperopic 

ablation (Fig. 17A) and oblate SA is induced after 
myopic ablation (Fig. 17B) . 

The epithelial smoothing model predicts that the 
minus cylinder ablation pattern of the Nidek EC-5000 

10 (Fig. 18A) achieves less correction with increasing 

smoothing action (Fig. 18B) . The clinical value of 
0.74 for the correction/ablation ratio of 0.74 is 
matched at s = 0.55 mm. The smoothing model also 
predicts an increased spherical equivalent (SE) change 

15 to 0 . 59 D at s = 0 . 55 mm. The simulated epithelial 

thickness change is shown on Fig. 18C. The deviation 
from intended correction is explained by marked 
epithelial thickening on the flat meridian at the outer 
edge of the OZ . There is also epithelial thickening in 

20 the center and thinning at the outside edge of the TZ. 

A believed understanding is that during corneal 
surface smoothing, epithelium thins over bumps or 
islands and thickens to fill divots or relative 
depressions. The mathematical model clarifies and more 
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precisely defines this belief. The model differential 
equation (Equation 14) states that epithelium thins 
over surfaces that are more convex and thickens over 
surfaces that are less convex. Convexity is quantified 
5 by the local Laplacian of surface height -V 2 h' . The 

local Laplacian is the sum of the curvature along two 
orthogonal tangential axes. It is also equal to twice 
the inverse of the radius of curvature of the best-fit 
sphere to the local surface. Fourier analysis shows 

10 that the differential equation is equivalent to a low- 

pass filter. The low-pass filter removes high- 
frequency (small-scale) undulations on the surface. 
The mathematical models agree with the believed 
understanding of epithelial smoothing. 

15 The model, in accordance with the present 

invention, postulates that epithelium migrates towards 
lower areas on the cornea, with height being defined 
locally along the line perpendicular to the surface. 
It is hypothesized that epithelial cell migration is 

20 inhibited by contact with neighboring cells, and net 

movement toward the lower side occurs because there is 
less contact on that side. Contact inhibition is a 
well known phenomenon and has biochemical basis in cell 
surface molecules. 
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Although the model is derived from postulated 
epithelial behavior, it can also describe the smoothing 
action arising from the surface tension of the tear 
film. In LASIK, the "draping" effect of the flap is a 
5 potential smoothing mechanism. In PRK, subepithelial 

"deposition of extracellular matrix material is another 
mechanisms for surface smoothing. In all of these 
mechanisms, less convex areas of the cornea are filled 
in. 

10 The model explains many of the observed phenomenon 

after LASIK and PRK. For example, the model, in 
accordance' with the present invention, predicts that, 
after myopic ablation, epithelium thickens in the 
center. This has been observed after myopic PRK and 

15 LASIK. The model, in accordance with the present 

invention, also predicts the epithelial hyperoplasia 
over the peripheral ablation area that occurs after 
hyperopic LASIK. 

The model, in accordance with the present 

20 invention, also explains regression or loss of 

treatment effect that has been clinically observe and 
linked to epithelial changes. The smoothing model 
correctly predicts more regression should occur after 
hyperopic treatment than myopic treatment. 



10 
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The minus cylinder ablation pattern on the Nidek 
EC-5000 (Fig. 18A) produces more spherical shift and 
less astigmatism correction than expected. This is 
correctly predicted by the epithelial smoothing model. 
5 A quantitative match with the correction/ablation ratio 

is achieved at s = 0.55 mm. This agrees well with the 
0.63 ,mm value for s estimated from the hyperopic 
simulations . 

The model, in accordance with the present 
invention, also correctly predicts that oblate 
spherical aberration (SA) would be induced by myopic 
ablation . 

An average value for s of 0 . 5 mm was determined. 
By averaging the estimates based on myopic, hyperopic, 
15 and cylinder ablations, the effects of extraneous 

variations may be reduced. The s parameter fully 
characterizes the corneal surface smoothing model. The 
ablation map a(x,y) that pre-compensates for expected 
surface smoothing can be calculated using a 
20 deconvolution operation. 

A' - AH' /F' 

Equation 53 
"Deconvolution" is the reverse of the 
"convolution" operation that describes our smoothing 
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model. Surface smoothing is described by a 
multiplication (Equation 48A) of the ablation map with 
a low-pass filter F' in the frequency domain. 
Multiplication in the frequency domain is equivalent to 
5 a convolution operation in the spatial domain. 

Equation 53 describes a de-filtering operation in the 
frequency domain that removes the effect of smoothing 
filter F' to achieve the desired surface change AH' . 
Direct solution of Equation 53 in the frequency domain 
10 is often unstable because F' become zero at higher 

frequencies (division by zero is undefined and division 
by a small number amplifies noise) . Thus, Equation 53 
is generally solved using iterative deconvolution 
procedures. 

15 The epithelial smoothing model, in accordance with 

the present invention, suggests that the target change 
map Ah(x,y) should contain a transition zone (TZ) 
around the OZ where there is a gradual change in 
curvature and slope. An abrupt change in slope would 

20 lead to an undesirable spike in radial' curvature and 

the convexity term V 2 Ah' (x,y). Even a sharp change in 
curvature would lead to a step transition in the 
ablation. According to the model (Equation 48), it is 
important to design the transition zone so that the 
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ablation does not contain much spatial frequency above 
1/s (radians/mm) . The hyperopia treatment would 
require a wider transition zone because there will be 
more phases of transition as well as greater slope 
5 change compared to myopia correction. 

The ablation process, in accordance with the 
present invention, compensates for the effects of 
epithelial smoothing and thereby improve the accuracy 
of refractive correction and reduce undesirable 

10 secondary aberration. A healing-adjusted ablation 

design minimizes induced aberrations in laser 
refractive surgery. Ablation designs in accordance 
with the present invention can be used to correct the 
full range of refractive errors such as myopia 

15 (nearsightedness), hyperopia (farsightedness), 

astigmatism, and higher order aberrations. The present 
invention makes wavef ront-guided treatment of 
aberrations more effective by anticipating and 
correcting secondary aberrations. It enables accurate 

20 results at the primary treatment, which is much 

preferred by patients over the necessity of undergoing 
a second surgery. 



